In any imaging system, from desktop scanners to remote sensing systems, the basic process is the same: signals from the object mixed with noises are quantized to form the image of the object. An imaging system may be conceptually viewed as a spatial convolution filter that acts on the true signals of objects and outputs the altered signals stored in an image. This spatial filter is the system Point Spread Function (PSF), which is the combined result of the system electro-optical characteristics and scanner motions. If the noise terms are assumed to be unrelated to the system and the object, then the process of remote sensing can be characterized by a so-called image equation:I=PSF{circle around (x)}O+N Where {circle around (x)} is the convolution operator. In this rather simplified image equation, the term I is the quantity directly measured by the imaging system and recorded in an image, O is the quantity to measure and is related to I through the known deterministic PSF and the unknown random noise term N. The process of resolving O from the observation I is often called image restoration or image reconstruction.
If both object and image are represented in the same coordinate system (x, y), this object-image relationship can be further expressed by a spatial integral:
      g    ⁡          (              x        ,        y            )        =                    ∫                  -          ∞                ∞            ⁢                        ∫                      -            ∞                    ∞                ⁢                              h            ⁡                          (                                                x                  -                  ξ                                ,                                  y                  -                  η                                            )                                ⁢                      f            ⁡                          (                              ξ                ,                η                            )                                ⁢                                          ⁢                      ⅆ            ξ                    ⁢                                          ⁢                      ⅆ            η                                +          n      ⁡              (                  x          ,          y                )            where
f(x, y) is the spatial representation of the true object, or O=f(x, y)
g(x, y) is the spatial representation of the image, or I=g(x,y)
h(x, y) is the PSF of the imaging system
n(x, y) is the noise function of the imaging system
As stated above, the task of restoration is to “invert” the imaging process and restore the original object f(x, y), based on the measurement g(x, y) and knowledge about the system PSF and the noise term. However, f(x, y) is not known and could never be completely restored. Therefore, the practical task has been shifted to finding an estimate {circumflex over (f)}(x,y) for f(x, y), such that it satisfies certain restrictions and optimal criteria.
Restoration of the observed image to the original object truth can be handled by the increasingly popular spatial convolution approaches. The centerpiece of the spatial convolution technique is a spatial filtering kernel, referred to hereafter as a restoration kernel, r(x, y); the estimation {circumflex over (f)}(x,y) of f(x, y) is calculated through the spatial convolution of the image g(x, y) and the kernel r(x, y):{circumflex over (f)}(x,y)=g(x,y){circle around (x)}r(x,y)
The restoration kernel r(x, y) is often required to satisfy certain restrictions, such as non negativity (r(x,y)≧0), and/or normality (∫∫r(x,y)=1). As the restoration kernel is key to the restoration of an image, there is a need in the art for building a restoration kernel that is particularly suited to specific applications.
The Fourier transform is a powerful tool for handling equations involving convolution and correlation. By applying the Fourier transform on both sides of the image equation, and using the convolution property of the Fourier transform, the following is provided:F[g(x,y)]=F[h(x,y)]·F[f(x,y)]+F[n(x,y)]or, the equivalent image equation in the Fourier frequency domain:G(u,v)=H(u,v)·F(u,v)+N(u,v)where G(u,v), H(u,v), F(u,v) and N(u,v) are the Fourier transforms of g(x,y), h(x,y), f(x,y) and n(x,y). In principle, both image equations in the spatial and the frequency domains are equivalent, since the solution for one equation can be derived from the solution for the other, either through Fourier transform or Inverse Fourier transform.
Likewise, the spatial convolution restoration equation:{circumflex over (f)}(x,y)=g(x,y){circle around (x)}r(x,y)can be written in Fourier frequency domain:{circumflex over (F)}(u,v)=G(u,v)·R(u,v)where {circumflex over (F)}(u,v), G(u,v) and R(u,v) are the Fourier transforms of {circumflex over (f)}(x,y), g(x,y) and r(x,y) respectively.
In consideration that H(u,v) may be zero for some frequencies, and N(u,v) can be significantly different from zero, the conventional Wiener approach looks for a filter function w(x,y), and an estimate function {circumflex over (f)}(x,y) such that:{circumflex over (f)}(x,y)=g(x,y){circle around (x)}w(x,y)and
      E    ⁢          ⌊                        (                                    f              ⁡                              (                                  x                  ,                  y                                )                                      -                                          f                ⋒                            ⁡                              (                                  x                  ,                  y                                )                                              )                2            ⌋        =                    min        γ            ⁢              E        ⁡                  [                                    (                                                f                  ⁡                                      (                                          x                      ,                      y                                        )                                                  -                                                      g                    ⁡                                          (                                              x                        ,                        y                                            )                                                        ⊗                                                            r                      s                                        ⁡                                          (                                              x                        ,                        y                                            )                                                                                  )                        2                    ]                      =          ɛ      2      where {rs(x,y)} is the set of all possible linear stationary restoration filters. In this expression, the expectation E should be understood as an average over all instances of the random noise field n(x,y) and over the whole (x,y) space where f(x,y) has definition.
Let W(u,v) be the Fourier transform of w(x,y), and assuming that the scene (object), the image, noise and filters are stationary, it can be proven that:
            W      ⁡              (                  u          ,          v                )              =                            R          Wiener                ⁡                  (                      u            ,            v                    )                    =                                    H            *                    ⁡                      (                          u              ,              v                        )                                                          H            ⁡                          (                              u                ,                v                            )                                ⁢                                                  2                        ⁢                                          +                                                      S                    n                                    ⁡                                      (                                          u                      ,                      v                                        )                                                              /                                                S                  f                                ⁡                                  (                                      u                    ,                    v                                    )                                                                          or            W      ⁡              (                  u          ,          v                )              =                  1                  H          ⁡                      (                          u              ,              v                        )                              ⁢                                              H            ⁡                          (                              u                ,                v                            )                                ⁢                                  2                                                          H            ⁡                          (                              u                ,                v                            )                                ⁢                                                  2                        ⁢                                          +                                                      S                    n                                    ⁡                                      (                                          u                      ,                      v                                        )                                                              /                                                S                  f                                ⁡                                  (                                      u                    ,                    v                                    )                                                                        where Sn(u,v)=Snn(u,v) is the noise power spectrum, and Sf(u,v)=Sff(u,v) is the scene (object) power spectrum. The ratio Sn(u,v)/Sf(u,v) can be viewed as a system noise-signal ratio. It can be shown that the modulus and phase of the Wiener filter are such that the Wiener filter is an inverse phase filter and a modified direct inverse modulus filter. It effectively controls the noise enhancement of the inverse filter. Additionally, the Wiener filter is a bounded filter. Because of this property, the Wiener restoration kernel has become more acceptable in image processing.
One prior art approach is by L. Wood for Landsat remote sensing data. In her 1986 dissertation entitled “Restoration for Sampled Imaging System,” incorporated herein by reference, Wood proposed a Wiener kernel construction approach for restoring Landsat Multispectral Scanner (MSS) data. Wood's approach has a number of weaknesses.
(1) A key assumption in Wood's approach is total separability. That is, every single component in the expression of the Wiener kernel is separable for x and y. However, even if all Wiener filter components are separable, it does not necessarily follow that W(u,v) is separable. To overcome this difficulty, the approach assumes that the Wiener filter itself is separable. This total separability assumption results in a suboptimal solution because a kernel defined this way is only an approximation of the true Wiener solution.
(2) Wood requires a number of estimations and assumptions to be made that may not hold true for all applications. For example, Wood's approach replaces the object power spectrum Sf(u,v) using the image power spectrum Sg(u,v). This manipulation may create some problems, because this manipulation implicitly assumes that the system noise-to-scene (or object) ratio is: 1−H2(u,v), which may not always hold true
(3) Wood's approach is strongly model dependent. For instance, the forms of the Optical Transfer Functions (“OTFs”) of the optical, detector, and the electro-filtering systems were assumed to satisfy respective separable mathematical models. The autocorrelation function of the object was assumed to satisfy an exponential function; the noise power spectrum of the noise was assumed to be one. The optimal nature of the derived filter can degenerate because so many specific function forms were hypothesized, especially the auto-correlation function of ground truth.
(4) Wood's approach relied on a particular scene and, therefore, it is not quite adaptive.
Thus, it can be seen that there is a need in the art for a new adaptive Wiener filtering algorithm for image data.